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Orbital measurements of neutrons by the Lunar Exploring Neutron Detector (LEND) onboard the Lunar 
Reconnaissance Orbiter are being used to quantify the spatial distribution of near surface hydrogen (H). 
Inferred H concentration maps have low signal-to-noise (SN) and image restoration (IR) techniques are 
being studied to enhance results. A single-blind, two-phase study is described in which four teams of 
researchers independently developed image restoration techniques optimized for LEND data. Synthetic 
lunar epithermal neutron emission maps were derived from LEND simulations. These data were used as 
ground truth to determine the relative quantitative performance of the IR methods vs. a default 
denoising (smoothing) technique. We review and used factors influencing orbital remote sensing of 
neutrons emitted from the lunar surface to develop a database of synthetic “true" maps for 
performance evaluation. A prior independent training phase was implemented for each technique to 
assure methods were optimized before the blind trial. Method performance was determined using 
several regional root-mean-square error metrics specific to epithermal signals of interest. Results 
indicate unbiased IR methods realize only small signal gains in most of the tested metrics. This suggests 
other physically based modeling assumptions are required to produce appreciable signal gains in 
similar low SN IR applications. 

Published by Elsevier Ltd. 


1. Introduction 

NASA’s Exploration Systems Mission Directorate developed the 
now orbiting Lunar Reconnaissance Orbiter (LRO) to quantify 
resources for future human activities on the lunar surface (Chin 
et al., 2007). The operating Lunar Exploration Neutron Detector 
(LEND), onboard LRO, was included to provide direct geochemical 
evidence for water on the Moon (Mitrofanov et al., 2008). 
However, the efficiency in which neutrons can be detected in 
lunar orbit is limited by the coupled effects of instrument design 
constraints and the low lunar neutron fluences at orbital altitudes 
(Feldman et al., 1993). This results in higher uncertainties or 
lower signal-to-noise (SN) at the detector which is offset by 
accumulating samples over a surface region using multiple orbits. 

Low SN in geochemical maps has been addressed in similar 
analytical efforts on data from the precursor Lunar Prospector 
(1998) mission's neutron spectrometer (LPNS) and gamma-ray 
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spectrometer (LPGRS) experiments that used image restoration 
(IR) techniques to deblur and denoise maps (Elphic et al., 2005; 
Lawrence et al., 2006; Eke, 2001). In this paper we evaluate the 
potential performance of several IR methods for unbiased 
restoration of expected lunar neutron emission derived from 
LEND (Mitrofanov et al., 2008). This includes the results of a 
single-blind study of IR method performance via direct compar- 
ison to synthetic neutron maps and against a default denoising 
technique designed to quantify signal gains specific to IR 
deblurring. 

While LEND is similar to LPNS in its operations and scientific 
objectives, there are significant differences in instrument design, 
expected mapping resolutions and SN that warrant independent 
evaluation of IR methods. For our single-blind study, four out of 
twelve invited teams of IR researchers accepted collaboration 
invitations. Teams configured and submit a total of four IR 
methods plus a default Gaussian denoising (smoothing) method 
for evaluation. The blind study was performed in a two-phase 
process similar to Nesbitt (2004). A preliminary training phase 
was used to configure and optimize IR methods which provided 
both “true" and degraded maps. In the subsequent blind study 
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phase only “true" maps were provided for restoration {Ivatury 
and McClanahan, 2009}* The primary research objectives are to 
identify optimal 1R methods for LEND map restoration as well as 
to quantify the expected improvement vs* the default map 
denoising method* The IR methods and teams include: 

1) Iterative Gaussian Smoothing (GSFC-IGS), T, McClanahan, 

NASA/Goddard Space Flight Center. 

2) Regularized Deconvolution (GSFC-UM-RD) r V. Ivatury, Univ* of 

Michigan/GSFC 

3) Weiner (UMD-Fourier), D. Usikov, G. Nandikotkur, G. Milikh, 

R. Sagdeev, Univ, Maryland (MD). 

4} Conjugate Gradient (UMD-CG), Usikov, G* Nandikotkur, 

G, Milikh, R. Sagdeev, Univ* MD. 

5) Pixon (PIXON), Pixon LLC*, R. Peutter. 

The duration of the paper reviews the necessary background in 
orbital lunar neutron detection and LEND specific instrument 
configuration used to formulate a model of the expected 
epithermal maps* A short description of the IR methods that 
were selected is provided as well as a description of the protocols 
that were used for the two-phase blind study. This review also 
includes the design and formulation of synthetic LEND maps, 
training, parametric configuration of methods and IR evaluation 
metrics* The results section quantifies IR method performances and 
signal gains against the default denoising (Gaussian smoothing) 
technique, 

2. Background 

The LRO mission was created by NASA in response to the 
Presidents 2004 Vision for Space Exploration and the tantalizing 
evidence of hydrogen deposits found at the lunar poles by LPNS 
(Arnold, 1979; Feldman et ah, 1999b; Feldman et aL, 2000). LRO 
successfully entered lunar orbit in June 2009 with scientific 
objectives that include locating and quantifying volatile deposits, 
e.g, water, and other potential resources that would support 
future human exploration of the surface. LEND is measuring lunar 
surface neutron emissions to refine the spatial distributions and 
quantity estimates of the putative polar H deposits identified by 
LPNS (Mitrofanov et aL, 2007}* 

Hydrogen on the Moon has been postulated for over 30 years 
(Arnold, 1979), Though originally thought to be dry, recent re- 
analysis of Apollo-15 lunar basalts collected at mid-latitudes 
indicates ambient H concentration may approach 745 ppm-0*07% 
(Saal et al., 2008). However, LPNS low resolution evaluations of 
the polar region indicated H concentrations may approach —3-5% 
and that enhanced concentrations may exist in smaller spots at 
the poles (Feldman et al*, 1998a). It is inferred from these results 
that H enhanced regions are accumulated through the combined 
effects of a regional H budget driven by variably offsetting factors 
governing depositional and entrainment processes* H depositional 
processes include periodic cometary and meteoritic bombard- 
ment and/or isotropic influx of solar wind protons (Vondrak and 
Crider, 2003). 

H accumulations can be regionally influenced by various 
geophysical entrainment effects that constrain the rate of H 
sublimation, At extremes some regional sublimation processes 
may be completely attenuated in polar topographic depressions, 
where persistent low polar inclination —1*5' to the sun produces 
low and permanent shadowing conditions (PSR)’s with persis- 
tently low temperatures < 1 00" K (Neubert et al, 2005; Vasavada 
et al., 1999). These conditions, over time are postulated to 
minimize H losses from sublimation and assuming isotropic H 
influx, may lead to regional H accumulations* As a result, PSR's 


have been defined as the primary LEND surface targets (Mitro- 
fanov et ai* p 2010; McClanahan et al., 2009). 

Importantly, the interpretation of orbital neutron measure- 
ments only indicates the presence of H, not its geochemical form 
e.g. water (H 2 0), hydroxy] ((OH)^ 1 ) and protons (H + ). Neutrons 
are produced in the lunar regolith via isotropic influx of high 
energy galactic and solar cosmic rays. In this process, subsequent 
nuclear interactions induce a number of reactions including 
spallation of neutrons and y-ray’s in the top few centimeters of 
regolith. Neutrons propagate through the regolith interacting 
with atoms, thereby thermalizing (losing energy) until reaching 
thermal equilibrium and finally being absorbed in the regolith, or 
importantly, lost into space as a leakage flux. Rates at which 
neutrons thermalize (lose energy) in materials are dependent on 
the material density and atomic cross section of atoms encoun- 
tered during this process. Compared to most elemental constitu- 
ents of planetary materials, H has a high cross section to neutrons. 
As a result, the energy distribution of the neutron leakage flux 
reflects intensity deficits at epithermal energies commensurate 
with higher H concentrations (Feldman et aL, 1993), Lunar H 
concentrations are derived from defected neutron distributions 
cross referenced to radiation and neutron transport modeling that 
factors regolith elemental compositions with variable H, neutron 
energies and detector/geometric configurations, e.g. Monte Carlo 
N-Particle Transport (MCNPX) codes (Pelowitz, 2008)* 

The areal scale of most polar PSR targets is small <few km 
diameter and the LPNS uncollimated field of view (FOV) at low 
altitude is inferred to be larger ( — 44 km FWHM) than the 
diameter of many PSR's (Feldman et al., 1999a; Bussey et al., 
2003). From the areal disparity in the expected small area of PSR’s 
and the larger LPNS orbital FOV, it is inferred some localized 
hydrogen deposits may exceed LPNS estimates due to the 
effective averaging (blurring) of variable epithermal rate regions 
within the FOV. LEND's design enhances its spatial resolution 
using a passive neutron collimator containing neutron absorbing 
materials, 10 B and polyethylene, to effectively discriminate sur- 
face emissions of epithermal neutrons to within ± 5.6 C of the 
instrument boresite. From its nominal 50 km mapping altitude, 
this configuration yields an instrument FOV that subtends a 
— 10 km diameter surface region (footprint). LEND will primarily 
point nadir and integrate signals at 1 Hz rates using its 10 internal 
and external detectors. These detectors monitor the energy 
spectrum of low (thermal), medium (epithermal) and high (fast) 
energy neutrons emitted from the lunar surface and from the 
spacecraft* To increase SN, four identical co-aligned nadir pointing 
epithermal neutron detectors are implemented in the collimator, 
as well as a total of six detectors positioned internal and external 
to the detector system for tracking background, instrument and 
spacecraft induced neutron fluxes (Mitrofanov et al., 2008, 2009)* 

2*1. Image and epithermal neutron map reconstruction 

Image formation is a stochastic process that convolves a true 
signal f, degraded by instrumental or environmental blurring w, 
with additive noise functions n, to produce an output image g, 
where g~f®w + n (Peutter et al., 2005; Gonzalez and Woods, 
1993). Instrumental blurring processes degrade the original image 
f as a convolution cg> with the instrument PSF, w and statistical 
counting noise is added, n (Poisson). Image restoration techniques 
aim to reverse the blurring process w 1 and attenuate the noise n, 
to obtain an estimate of the original image, g In high SN 
applications many types of image restoration method have been 
shown to yield significantly improved results (Peutter et al., 
2005). However, image restoration becomes increasingly proble- 
matic with lower SN applications due to the combined effects of 
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statistical noise being variably spatially correlated, In these 
conditions spatial frequencies may intersect with real image 
frequencies and increase the difficulty in estimating an under- 
lying image formation model. Additionally, some restoration 
processes may amplify noise signals, yielding artifacts e.g. 
fluctuations (mottling), which strongly resemble real signals 
(McClanahan ef al., 2007; Ugbeme et al., 2007). 

Accurate map formation models for orbital neutron remote 
sensing are also difficult to estimate due to the low signal rates 
and variances in environmental and instrumental sampling 
conditions inherent in long duration mapping. Efforts to mitigate 
these issues and to optimize reconstructions include techniques 
to extend map formation models by use of other physically based 
constraints. However, the implementation of physical constraints 
in the 1R process is performed through either parametric biases or 
to establish prior physical limits for some pixel values (Elphic 
et al., 2007; Eke et al., 2009). Eke et al. claim signal gains in 
restoring LPNS data only with the implementation of physically 
based constraints, e.g. shadowing. Concerns with this approach 
are that the lunar polar hydrogen budget is presently poorly 
understood. Implementing these assumptions as restoration 
priors may neglect either the impact or appropriate weight of 
other geophysical processes. Results of non-biased IR techniques 
on LPNS data were equivalent to the same data smoothed with a 
Gaussian operator, underscoring the importance of accurate 
priors and understanding the ramifications of potential bias 
errors in low SN IR applications (No Free Lunch) (Wolpert and 
Macready, 1997), 

Fig. 1 illustrates LEND's neutron map formation process 
including required functions and parameters to configure the 
described IR methods. This includes surface emission and instru- 
mental processes e.g. PSF w, background estimation and techni- 
ques used to map neutron detections from orbit. Poisson counting 
statistics, n describe the arrival process of neutrons accumulated 
by the detector over a given region and defines the expected pixel 
uncertainty as sqrt(a). A consequence of this long term multi- 
sample mapping strategy for IR e.g. LPNS, is that as opposed to 
traditional higher SN IR models, each detector measurement 
accumulated to the map is uniquely a function of time, instrument 
operations, position and background. As a result, new variables are 
introduced to the map model e.g. altitude, solar flux, position, 
pointing which impact estimates of the PSF and statistical 
uncertainties needed for image reconstruction. For instance, 
altitude variance in combined measurements at a given pixel i, 


LEND collimated 3 He sensor 



Fig. 1. A single LEND 3 He collimated detector configuration and orbital factors 
defining a nadir pointing collimation process, lunar surface flux of neutrons 
incident the collimator exterior are passively discriminated using neutron 
absorbing materials. Collimation yields a pixel size — lGkm from 50 km 
altitude (E), 


adversely impacts both the instrument blurring function w and 
statistical uncertainty n estimates at f. As a result, accumulation of 
neutrons from variable sampling conditions to the map incurs 
additional systematic errors. 


3, Methods 

The IR performance study was implemented as a two-phase 
process using four independent research groups. Section 1. GSFC-1GS 
is defined as the benchmark "denoising only” approach and 
establishes the null hypothesis. Each IR method uses different 
numbers and types of parameters which must be configured for a 
given application domain. To configure and optimize these methods 
for LEND map restorations, a preliminary set of 12 synthetic LEND 
epithermal maps were provided, simulating north, south, polar (NP, 
SP) and randomized PSR spatial and intensity distributions. Training 
phase data was distributed to teams and included both degraded 
synthetic LEND maps g and ground truth maps f. Maps for both the 
training and blind study phases were produced using LRO mission 
planning ephemeris for the 382 day primary exploration mission 
with the assumption of continuous nadir pointing of the LEND 
instrument and 100% duty cycle. This virtual flight process 
effectively flew the LEND instrument mission over quasi-randomly 
derived lunar polar neutron emission models ( ± ) 80-90* factoring, 
ephemeris, 1 second integrations and LEND’s expected PSF 
(Mitrofanov et al„ 2009). 

A consequence of LEND's collimated design is that the exact 
surface emission points of detected neutrons are only statistically 
defined within the FOV. This spatial uncertainty constitutes the 
present basis of LEND f s symmetric, 2-D normal PSF, w, N{^= 0, 
t7 XJ ,= 1.67 km) from 50 km altitude. Where, w defines the spatial 
probability distribution of subtended pixels. The PSF is defined 
theoretically via the 5.6° aperture and collimator length, assuming 
complete discrimination of incoming epithermal neutrons incident 
to the external side of the collimator. Only epithermal neutron 
detection processes are assumed with statistical counting uncer- 
tainties and no degradation due to motion blur. The map formation 
process for training and blind studies is illustrated in the flow 
diagram, Fig, 2. Examples of the LEND map model components are 
illustrated in Figs, 1 and 3, where the LEND 382 day coverage 



Fig. 2. Flow diagram for generating syntherie lunar epithermal maps. 
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distribution over the NP region is depicted in Fig, 3a. A synthetic 
true lunar epithermal map with craters as depressions in the 
epithermal count rate (dark spots) f r is defined in Fig, 3b, Craters 
delimit sources of enhanced (H) that have lower epithermal flux 
rates than surrounding background. These regions deviate ran- 
domly up to 18% below the expected 0,88 cps collimated rate 
which is equivalent to the 4.5% suppression seen by LPNS at the 
lunar north pole. (Feldman et al., 1 998b). Crater locations, sizes and 
intensities were defined in (40) maps randomly, and in known (40) 
south polar and (40) north polar craters defined by International 
Astrophysical Union (IAU) conventions (Anderson and Whitaker, 
1982; Bussey and Spudis, 2004). 

Blurring of the image in Fig. 2 is performed and the product of 
the blurred time and blurred true image (t® w)*(f® w) provides 
the counts map upon which random Poisson statistical uncer- 
tainties n are added. Factors for randomly deriving the synthetic 
neutron epithermal maps used in training and single-blind study 
are listed in Table 1 . 


3,1. Map restoration methods: iterative Gaussian smoothing 
(GSFC-IGS) 

Gaussian smoothing (GS) has significant precedence in 
signal denoising applications and has been used to smooth 
published planetary neutron and gamma-ray geochemistry maps 


(Mitrofanov et al., 2002; Boynton et al., 2004). Iterative Gaussian 
Smoothing (IGS) has been shown to be an equivalent process to a 
diffusion process {Koenderink, 1984; Perona and Malik, 1987). 
This denoising method systematically smooths data by incremem 
tally attenuating the power in high frequency signals first and 
progressing towards attenuating the power in lower spatial 
frequencies Costa and Cesar (2001). In the GSFC-IGS approach, 
the detected signal, g is iteratively smoothed j times. For each 
iteration the map g, is convolved with a fixed size, unit volume 


Table 1 

Synthetic neutron emission map generation factors for training and blind study. 


Map size 

Map evaluation area 

Predefined North Craters (40) 
Predefined South Craters (40) 
Predefined Random Craters (40) 
Crater intensities 

Accumulated signal scale/382 day 
baseline 

Background intensity (blobby. variable 
smooth) 

Point spread function (field of view) 

Ephemeris 382 days 

Spacecraft poiar altitude variable 


1620 *620] pixels, pixels 1 km 2 
+ 80-90°. pole pt [310. 310) +303 km 
radius 

85 (size, position fixed per crater tD) 
129 (size, position fixed per crater ID) 
Random Size . Random Position 
{3,6, 9. 12. 15, 18] % decrease 
(random/crater) 
j 1.0, 1 75, 2.5j x Accum, cover, 
(random/crater) 

Mean— 0.88 cps. range 15% of mean 

LRO, Beckman-V4, 2007 
[30-60] km, mean-50 km 


a 


c 


Time Spatial Distribution, t 
Pole 

t = time 




Blurred, Noisy 
Count Rates, cps 


With 
Poisson 
Noise, g 


True f, r 
-l- Poisson, g 


3.43 




Blurred and Noisy, g 



Deconvolved Map, F 

True, Deconvolved 



Fig. 3. (a) (Top left) LEND lunar pole accumulated spatial coverage (t^time) distribution 80-90 7 for 382 days ephemeris using 1 km 2 pixels. Simulation projects the LEND 
FOV in 1 $ integrations to a base map. Double lobes in the distribution are due to lunar gravitational anomalies impacting LRO trajectories near the pole (P). (b) (Bottom left) 
Example true count race map f/t, time normalized with randomized placement and crater intensity (dark spots, variable intensity) blobby background, p-0.S8 cps. range 
< ±7.5%. Range-[0.57. 0.93] cps. Craters represent randomly defined, discrete decreases in epithermal count rates. [3-lSj%. Area >80' : latitude evaluated by 
methods =28841 3 pixels/map. (c) (Top right ) Example time normalized blurred and noisy map g/t , after convolution of true map (b) with LEND FOV and added Poisson 
statistical uncertainties from map used as input to IR methods. Range -{0.. 3.43) cps. Metrics compare reconstructed maps against the true map (b). (d) (Bottom right) 
Example of LEND epithermal map after image restoration (PIXON), f/t. Range -{0.57, 0.93), 
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2-D Gaussian smoothing kernel k, gj+i=gj®k, For each true 
training map f, the root-mean-square (RMS) error is quantified at 
each iteration j comparing true f and smoothed gj. For each 
training map t, the minimal RMS error is determined after each 
iteration. An average of the IGS RMS error training curves yields 
[11 iterations] as the optimal iteration setting; to denoise the 
blind study dataset (GSFC-ICS), 

3.1.1. Regularized IR: (GSFC-UM-RD) 

The GSFC-UM-RD method optimizes a regularized deconvolu- 
tion function (DECONVREG) available via the MATLAB® signal 
processing package. A regularization strategy minimizes the 
degree of noise amplification incurred from division of smalls 
0, high frequency Fourier coefficients common in Weiner type IR 
methods (Weiner, 1942), In this approach, inverse filtering is 
performed using a variant of the least squares method, 
rig-~w<g>f] 2 +a[pts>f] 2 . Regularization term p defines a discrete 
Laplacian distribution convolved with f yielding a smoothed 
function scaled by a (LRANGE), defining a spatially specific noise 
minimization constraint A second parameter in the function is 
known as the additive noise power, np of the image, np is 
calculated using a 2-D Fast Fourier Transform (FFT) of g, A center 
inverse Fourier transform is then applied to the blurred image g. 
Optimal settings for NP, LRANGE parameters were defined via 
gridded search of the discretized DECONVREG input parameter 
space. Optimal input parameter settings |32, 993], 

3.1.2 ' Weiner (UMD-Fourier) 

A fast Fourier transform (FFT) based method, developed by 
{Press et aL, 2008 ; Weiner, 1 942) is well known for speed, and has 
application to signals corrupted by noise and coordinate-invariant 
point spread function w (Section 3.1), The simplicity of the 
method lies in the fact that the FFT of the detected image is the 
product of FFTs of the blurring function w and true image f, 
FFT(g) — FFF(w)*FFT(f) thereby enabling inversion FFT(f) = FFT(g)/ 
FFF(w), yielding an estimate of the true image (FFT"' 1 ), f. 
However, the inversion becomes problematic in the presence of 
noise (low SN), and small high frequency coefficients in the 
denominator that can amplify noise. We address the issue of small 
FFT coefficient division using additional frequency dependent 
terms in the denominator [a, b\ (a variety of the Wiener filter) as 
FFT{g)/|FFT(w)+ci+b*FFr(f')]. f' is an estimate of the true image. 
This approach optimizes model parameters using the trial images 
via gridded search of the input parameter space and obtained the 
following settings [1e-04, le-05] (Rosenfeld and Kak, 1982). 

3.1.3, Conjugate Gradient (UMD-CG) 

Performance of the UMD-Fourier method is dependent on 
event statistics in the individual pixels. These problems are 
treated using a Conjugate-Gradient method (CG) (Nguyen et al., 
2001). Reconstruction is performed by minimizing the least- 

square-method (LSM) sum, S= j(w T f-g) T £~ 1 (w T f“g) which 
translates to solving the matrix equation Af-w T g, where the 
A^w t ^ _1 w, and terms w, J2* f and £ are the point spread 
function, the covariance matrix of event statistics (in our case — a 
diagonal one), the true(unknown) and convolved images in 
vectorized form. Although the solution is a simple matrix 
inversion, its computational complexity, O(620 2 ) 2 precludes 
processing via conventional computer hardware. 

Alternatively, gradient descent methods such as CG can 
mitigate the inversion limitation by iteratively adjusting method 
parameters via discrete estimation of parametric gradients. 
Methods converge in a limited number of iterations assuming 
initial parametric estimates are proximal to the global solution. 
For low SN, a regularization term is included, e.g. A=rW T J2“ 1 w+ 


yU where y is a positive regularization coefficient for noise 
suppression, L is a Laplacian circular distribution. The Least 
Square Method (LSM) sum becomes S=^(w T f-g) T r mt (w r f- 
g) + yL£, Simulations with the trial images yielded optimal 
reconstructions using [600 iterations, y = 1000], 

3.1 A. Pixon LLC (PIXON) 

The Pixon* 1 method is a picture-information-based method 
that performs image reconstruction under the constraint that the 
model for the image should be the simplest model consistent with 
the data {Puetter et al., 1995, 2005). PIXON gains its power by 
imposing a minimum complexity constraint on the reconstructed 
image. Unlike many other competing techniques, the minimum 
complexity constraint strongly suppresses unwanted image 
artifacts, as such artifacts are unneeded in fitting the data, and 
hence are by definition removed. PIXON proceeds by localized 
smoothing of a minimized set of independent map sub-regions 
(Pixon elements). Formally, this is described as an integral over 
the map,I(y) = / dzK(y,z)<p(z), where I is a pseudo-image, K is a 
positive valued kernel function, <p is a spatial scale adaptive 
regularizing term sensitive to localized noise estimates and limits 
high frequency fluctuations at scales <K. PIXON regions are 
dynamically size adaptive as a function of local uncertainties. 
Convergence to a solution comprises a parallel optimization 
seeking largest kernel functions K and minimal least squares fit of 
the image data for each pixon. 

There are several options for implementation and configuring 
PIXON parameters. For this research INIT_ANNEAL and ANNEAL 
parameters were optimized in the training process using a 
gridded search of the 2-D parameter space [5e~6, 6e-4]. Maps 
were preprocessed prior to PIXON processing using a spatially 
adaptive flux preserving filter to identify regionally variable SN 
for an optimally smoothed result. 


4, Results: RMS error ±80-9(P 


The main goals of this research are to evaluate and quantify the 
denoising and debiurring properties of IR methods (Section 3.1) 
on expected LEND epithermal neutron maps. This includes 
(1) identify IR methods applicable to LEND (2) their parametric 
configuration (3) potential for induction of artifacts and 
(4) quantify the error reduction specific to IR vs. the default 
GSFC-IGS denoising approach. Metrics for method evaluation 
entail comparing the database of 120 restored maps f against the 
known synthetic maps f, by averaging the RMS error of map i 
regions/signals of interest reg, with pixels in regions j, Eq. (1). 
Secondary metrics include analysis of residual error distributions 
and pixels quantified as statistical outliers 


Regi'orc_RMSError = 


120 


E 


J - l 


( 1 ) 


Table 2 states each methods mean RMS error results averaged 
over the 120 blind study maps. Columns indicate region specific 
evaluations: (1) Global RMS error for ±80-90 y region followed by 
the decomposition of each map i into (2) Background, (3) Crater 
regions. Three best performing methods in each metric are color 
coded by rank, from lowest to highest RMS error [ 1 = green, 
2=ye//ow, 3=red]. Uncolored cells indicate rank positions >3rd. 
Top performing methods: (1 ) UMD-CG, (2) PIXON and {3) GSFC-IGS. 
However, PIXON performed nearly as well with first differences 
only noted at the 7th mean RMS decimal position. 

An average of 80% of all map area is allocated to background 
and 20% to craters, where each region is characterized as follows: 
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Table 2 

Method results: mean RMS error and method rank for 120 maps for 10° estimates 80-90’ and method rank in set (1 green (best), 2 yellow, 3 red) (column 1). 10’ results and 
rank using Minkowski L2 metric (Duda et al., 2001) (column 2). Background only mean RMS error (column 3), in crater regions mean RMS error (column 3). 


Method \ Metric 

10-degree 

MRMS 

10-degree 
L J Norm 

Background 

MRMS 

Crater 

MRMS 

Category 

Rank 

GSFC-IGS 

0.0211180 

11.3412 

0.0186679 

0.0286845 

1 

GSFC-UM-RD 

0.0213130 

11.4450 

0.0190138 

0.0284829 

2 

UMDCG 

0.0210508 

11.3052 

0.0178985 

0.0302985 

3 

UMD Fourier 

0.0217578 

11.6848 

0.0183534 

0.0316352 


PIXON 

0.0210573 

11.3086 

0.0175651 

0.0310757 



(2) Background: low amplitude spatial frequencies emulating the 
expected dynamic range and variation of lunar epithermal surface 
flux. Best performing methods identified (1) PIXON, (2) UMD-CG 
and (3) UMD-Fourier. Only small differences were noted between 
these results < 0.001 cnts/sec, RMS. In background regions the 
dominant spatial frequency is near the direct current (DC) 
coefficient. Amplitudes of other higher frequency coefficients 
are comparatively small so, convergence to the DC baseline is 
expected for both IR methods and GSFC-IGS. As such, GSFC-IGS 
performed nearly as well as comparative restoration methods, in 
4th place, with small first differences between GSFC-IGS and 
UMD-CG at the 4th decimal place [0.021 1 18, 0.0210508] cnts/sec. 

(3) Craters: Table 2, column three defines results of evaluation of 
localized small spatial scale depressions in the epithermal 
neutron flux. Cratered regions are of primary interest for NASA’s 
manned lunar exploration programs due to the potential for 
exploiting localized quantities of (H) as a resource. These regions 
are inferred to occur primarily in regions of low illumination, but 
are defined here as cratered regions with variable sized depres- 
sions in the lunar surface epithermal neutron flux with variably 
high spatial frequencies at crater edges. 

All method performances were worse than their (1) Global 
RMS errors, >0.007 RMS. Relative weighting of cratered and 
background results are consistent with the contributions factoring 
area-weights in the 10° latitude band results. Higher RMS errors 
are seen due to the overlap of higher spatial frequencies in crater 
rims with noise frequencies at lower latitudes and are considered 
more difficult spatial frequencies to reconstruct. The ranges of 
method results indicate performance was consistent across the 
study, range < 0.003 RMS. Best performers in this category were 
(1) GSFC-UM-RD, (2) GSFC-IGS and (3) UMD-CG. Importantly, 
with the exception of GSFC-UM-RD, results indicate that restora- 
tion methods did not improve significantly upon GSFC-IGS in 
cratered regions. Signal gains in cratered regions attributable to IR 
methods are of significant interest in this study due to the 
expected higher H concentrations and potentials for human 
exploration. Results in these higher spatial frequency regions 
did not identify significant gains vs. the IGS methodology. 
Additionally, UMD-CG and PIXON, use of a smoothing preproces- 
sing step prior to the restoration step may have incurred 
additional unrecoverable signal losses in these regions. 

4.1. RMS error vs. latitude 

Fig. 4 and Table 3 illustrate the RMS error results as a function 
of latitude. Map SN increases as a function of higher latitude via 
the combined effects of the time, LRO mission’s polar crossing 
trajectory, altitude and the expected overlapping of the LEND 
footprint over the poles. As a result the accumulated coverage is 
highest within ^3° of the pole where the greatest degree of 


RMS Error vs Latitude 



Fig. 4. Five IR methods latitude dependent mean RMS error averaged over 120 
blind study images. 


ground coverage overlap occurs. The accumulated coverage 
decreases significantly at latitudes below ~87°, thereby decreas- 
ing the SN, Fig. 3a. To evaluate the impact of variable SN on 
method restorations, the map is decomposed as a function of 20 
concentric latitude bands, 0.5 9 in size (15 pixels wide) in which 
average RMS error is calculated, Eq. (1). A consequence of this 
binning of results is that each latitude band area decreases 
towards the pole yielding unequal weighting of bands in Fig. 4. 

Moving from low to high latitude (outside map to center), 
80-82° indicates PIXON and UMD-CG methods have superior 
reconstruction performances. This region constitutes ~36% of the 
polar region evaluated. At 80°, images are dominated aerially by 
background regions with exceptions in some included craters in 
the randomized maps. As a result, due to lower spatial frequencies 
UMD-CG and PIXON performed best, which was consistent with 
both the overall and background results described in Table 2. As a 
fraction of area NP, SP craters are more prominent in the higher 
latitudes, which decreases the relative distribution of background 
to cratered pixels in these latitudes. Moving poleward to 81°, all 
performances are noted to improve, method performances 
intersect and RMS error decreases due to the higher SN 
approaching the poles. GSFC-UM-RD, GSFC-IGS results are 
equivalent, and alternate as the best performers. Trend disconti- 
nuities are indicated at 84.5° and 86.5°, which are correlated with 
the generally higher occurrences of fixed craters for NP, SP maps. 
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Table 3 

Mean RMS errors for 0.5 s latitude bands ± 80-90' indicating method performance and relative rank for restoring 120 blind study epithermal maps. Region from 80-82.5° 
constitutes ~50% of evaluated area. 


Latitude 

GSFC-IGS 

GSFC-UM- 

RD 

U-MD-CG 

UMD 

Fourier 

PIXON 

Band Area 
(fraction) 

Rank 

80 

0.0214532 

0.0215892 

0.0200488 

0.0208074 

0.0192914 

0.09741 

1 

80.5 

0.0216799 

0.0220767 

0.0210381 

0.0219247 

0.0205457 

0.09251 

2 

81 

0.0226211 

0.0230121 

0.0223008 

0.0230952 

0.0220945 

0.08753 


81.5 

0.0224759 

0.0229242 

0.0223110 

0.0230088 

0.02234 11 

0.08233 


82 

0.0223501 

0.0226315 

0.0222680 

0.0230125 

0.0223495 

0.07746 


82.5 

0.0221881 

0.0224350 

0.0222824 

0.0230964 

0.0224117 

0.07270 


83 

0.0214667 

0.0217276 

0.0215741 

0.0222858 

0.0217145 

0.06750 


83.5 

0.0213115 

0.0214573 

0.0214711 

0.0221484 

0.0216352 

0.06241 


84 

0.0208043 

0.0210435 

0.0210271 

0.0216268 

0.0212538 

0.05756 


84.5 

0.0210452 

0.0211569 

0.0215559 

0.0222956 

0.0218656 

0.05260 


85 

0.0203602 

0.0204726 

0.0208171 

0.0214978 

0.0211061 

0.04755 


85.5 

0.0192660 

0.0192408 

0.0197538 

0.0204231 

0.0200792 

0.04251 


86 

0.0181739 

0.0182990 

0.0186240 

0.0191865 

0.0189104 

0.03734 


86.5 

0.0191531 

0.0190952 

0.0199612 

0.0204816 

0.0203486 

0.03258 


87 

0.0177555 

0.0178467 

0.0184375 

0.0188290 

0.0187810 

0.02751 


87.5 

0.0160259 

0.0160105 

0.0166219 

0.0171157 

0.0168565 

0.02254 


88 

0.0151404 

0.0151743 

0.0158315 

0.0163040 

0.0160822 

0.01741 


88.5 

0.0148233 

0.0146576 

0.0157896 

0.0163187 

0.0160251 

0.01254 


89 

0.0145032 

0.0143719 

0.0156416 

0.0159537 

0.0158649 

0.00751 


89.5 

0.0158328 

0.0163019 

0.0172973 

0.0168458 

0.0177089 

0.00250 



Error trends are elaborated in Table 3 where, the top two 
performers in each latitude band are highlighted [l=green, 
2 = red]. 


4.2. Residuals 

Residual error distributions are derived by subtracting the true 
map from the restored map, r=gf-f. Histograms of r pixels 
provide the basis for metrics that evaluate the globalized IR 
reconstruction quality as well as the potential induction of map 
artifacts. Assuming residuals are normally distributed, a Gaussian 
fit of the r distribution provides both the mean [x r and variance <j 2 t 
terms. Binning of the residual errors was performed at 0.001 cps 
intervals. As residual pixels are identical in weight and total 
numbers, (~3.46 x 10 7 total pixels in the study), method residual 
error histograms have identical area. Direct comparisons of 
method residual distributions are performed as follows: (1) the 
residual error distribution mean fi r deviation from 0. quantifies 
the existence of bias in reconstructions. (2) Gaussian fit of the r 
distribution yields, amplitude and variance a 2 terms which are 
inversely related. Optimally, for best IR performers the fit 
amplitude is increased and variance is decreased. (3) Artifact 
induction is inferred by greater numbers of high error pixels. To 
quantify this error for each method, the integral error in r 
histogram tails (a) outside ± > 3cr r is derived. <r r from the 


(PIXON) method (lowest variance) was used as the baseline, 
°>ixon = 0-0163 cps. 

Fig. 5 illustrates fits of each methods residual error distribu- 
tion. Results of the bias study (1) indicated low bias for all 
methods < 0.0022-0.2% of the mean detector count rate 0.88 cps. 
Best performing methods were (1) GSFC-UM-RD, (2) GSFC-IGS 
and (3) UMD-CG. The amplitude/variance metrics are coupled due 
to their mutual dependence. These results identified: (1) PIXON, 
(2) UMD-CG and (3) GSFC-IGS as having the highest amplitude/ 
lowest fit variance. Values for each methods residual error metrics 
are stated in Table 4. 

Table 5 defines method performances in the analysis of 
residual tail areas. In each case the Tail Hi areas were 
approximately twice the areas of the Tail Lo areas. This is 
assumed due to Poisson statistics generating fewer 0. cut-offs 
for pixel counts. Relative areas in these regions are approximately 
Tail Lo: 1.3% and Tail Hi: 2.4% of the total pixels processed for all 
methods. For normally distributed data, each tail above ±3cr 
should contain approximately 0.15% of all pixels processed 
indicating that in all methods tail areas are above the expected 
area and not normally distributed. Each methods tail areas are 
integrated above and below the two tail thresholds established 
using the PIXON baseline positions + 3<7 PIX on- For Tail Lo) 
(1) UMD-CG, (2) UMD-Fourier and (3) PIXON had the lowest tail 
areas. For Tail Hi) GSFC-UMD-RD had the lowest area followed by 
GSFC-IGS and UMD-CG. A third metric combining tail areas was 
established after noting the larger numbers of high intensity pixel 
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Table 4 

Moments from Gaussian fit of residual error. Mean indicates existence of bias as deviation from 0. High amplitude and low sigma indicate higher number of pixels with low 
error. Rankings indicate all methods displayed low bias < 0.0023 cps, optimal GSFC-UM-RD. PIXON had highest amplitude and lowest variance followed by UMD-CG. 


Residual Error 
Analysis 

Method 

Amplitude 

Mean 

Sigma 

Category 

Rank 


GSFC-IGS 

748838 

-0.0015015 

0.0176412 

1 

GSFC-UM-RD 

735330 

-0.0009466 

0.0180208 

2 

UMDCG 

775114 

-0.0017858 

0.0169150 

3 

UMD Fourier 

738976 

-0.0022218 

0.0178038 


PIXON 

797804 

-0.0019891 

0.0163151 


Gaussian Fit of Method Residual Error Distributions 

900000 
800000 
700000 
600000 
| 500000 
'£ 400000 
300000 
200000 
100000 
0 

-0.080 -0.060 -0.040 -0.020 0.000 0.020 0.040 0.060 0.080 

Residual Error (cps) 

Fig. 5. Gaussian fits of method residual errors. Means indicate low bias for all methods. Variances indicate PIXON has lowest variance (also greatest peak amplitude), 
however examination of residuals in regions ± 3 a indicate IGS smoothing has the lowest tail area in the high error region. 



Table 5 

Integration of residual tail areas using (a) Tail Lo = residual errors < -0.048 cps, UMD-CG indicated best performance, (b) Tail Hi residual errors > 0.048 cps. GSFC-UM-RD 
indicated best performance. Combined tail areas indicated IGS performed best. 


Tail Error 
Analysis 

Method 

Tail Area 
Lo 

Tails Area 
Hi 

Total 

Areas 

Category 

Rank 


GSFC-IGS 

431880 

728209 

1160089 

1 

GSFC-UM-RD 

462744 

717249 

1179993 

2 

UMD CG 

393927 

802467 

1196394 

3 

UMD Fourier 

401568 

900394 

1301962 


PIXON 

405688 

844369 

1250057 


errors in the Tail Hi) metric. (3) Tail Totals yielded (1) GSFC-IGS, 
(2) GSFC-UM-RD and (3) UM-CG with the respectively lowest 
total tail area. 


5. Conclusions 

This paper described a two-phase blind study evaluation of 
four IR methods configured to evaluate and quantify restored 
lunar epithermal neutron flux maps. Background in orbital remote 
sensing of neutrons and LEND related factors was described with 
context in image restoration research. Comparison of IR methods 
was performed using a test database of synthetic LEND maps, 
derived from signals modeled from existing research and LEND 
instrument factors. IR methods were directly evaluated against a 


traditional smoothing process, GSFC-IGS to quantify the expected 
gains specific to IR deblurring vs. denoising. A 12-map training 
phase was included prior to the blind study to familiarize teams 
with the dynamics of LEND signals and facilitate method 
optimization. Teams optimized their IR methods on the training 
dataset via gridded search or simulated annealing methodology 
by evaluating method performance. GSFC-IGS was trained by 
establishing the optimal number of smoothing iterations (11) on 
the training set that minimized the RMS error. 

The blind study phase comprised each teams reconstructing 
the 120 degraded map database using their respective training 
configurations. 

Three types of IR method metrics were described (1) RMS error 
±80-90°, (2) RMS error vs. latitude and 3) residuals. In (1) only 
small gains were realized for UMD-CG and PIXON over 
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GSFC-IGS, < 0.0001 cps RMS. Relative improvements were 
approximately 2.4% of the expected LEND 0,88 cps collimated 
count rate. It is inferred this is primarily due to the larger relative 
distributions of background vs. crater area, where background 
regions are more homogenous in intensity. However, for cratered 
regions of interest, GSFOUM-RD, GSFC-IGS performed slightly 
better (0.02848, 0,02868] cps than UMD-CG, PIXON [0.0303, 
0,031 1 ]. (2) Map latitude analysis indicated method performance 
varied as a function of the relative distributions of background 
and cratered regions which varied with latitude. PIXON and UMD- 
CG performed best in map outer latitudes characterized by lower 
SN and higher relative numbers of background pixels, GSFC-UM- 
RD and GSFC-IGS performed better in map polar regions with 
higher SN and higher relative numbers of cratered pixels. Fig. 4, 
This result may alternately reflect differences in training method 
optimization. (3) Residual error analysis via Gaussian fitting of 
residuals identified PIXON with the best performance at mini- 
mizing fit variance / maximizing amplitude [0.016, 797804], 
However, the analysis of method tails indicates the GSFC-IGS 
method best minimized the existence of high error map pixels 
> +3 (7 vs. reconstructions, < 0,04% tail area difference. 

In summary, only small inconsistent signal gains in some 
metrics was realized by unbiased IR methods over the default 
GSFC-IGS, specifically in the global metrics. Importantly, in more 
detailed regional analysis (craters) with higher "true" spatial 
frequencies, three of four iR methods performance lagged GSFC- 
IGS in these critical signals. Similarly, all IR methods induced 
greater numbers of high error pixels (artifacts) vs, the default 
GSFC-IGS. Finally, IR methods were optimized with ground truth 
{generally not available in most applications) and yielded only 
small inconsistent gains vs. GSFC-IGS for the low LEND SN. These 
results would appear to limit the effectiveness of unbiased IR 
methods in similar low SN applications. However, methods 
incorporating other physically based constraints into low SN IR 
methods may have potential in improving these results but may 
incur the risk of additional reconstruction errors due to the 
incorporation of new image modeling assumptions. 

Future research will study the use and application of physically 
based modeling assumptions for LEND map reconstructions. 
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